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Abstract 

We review the existing mathematical models which describe physicochemical mecha- 
nisms capable of producing a symmetry-breaking transition to a state in which one chi- 
rality dominates the other. A new model is proposed, with the aim of elucidating the 
fundamental processes at work in the crystal grinding systems of Viedma [Phys Rev Lett 
J2> 94, 065504, (2005)] and Noorduin [J Am Chem Soc 130, 1158, (2008)]. We simplify the 

qt 1 model as far as possible to uncover the fundamental competitive process which causes the 

symmetry-breaking, and analyse other simplifications which might be expected to show 
symmetry-breaking. 
> 

^h ' 1 Introduction 

rn ; 

A significant stage in the formation of living systems was the transition from a symmetric 
chemistry involving mirror-symmetric and approximately equal numbers of left- and right-handed 
chiral species into a system involving just one-handedness of chiral molecules. 

In this paper we focus on mathematical models of one example of a physicochemical system 
which undergoes such a symmetry-breaking transition, namely the crystal grinding processes 



O 

o 



X 



investigated by Viedma [29] and Noorduin et al. [21] . which have been recently reviewed by 



McBride & Tully [18]. Our aim is to describe this process by way of a detailed microscopic 
model of the nucleation and growth processes and then to simplify the model, retaining only the 
bare essential mechanisms responsible for the symmetry-breaking bifurcation. 

We start by reviewing the processes which are already known to cause a symmetry-breaking 
bifurcation. By this we mean that a system which starts off in a racemic state (one in which both 
left-handed and right-handed structures occur with approximately equal frequencies) and, as the 
system evolves, the two handednesses grow differently, so that at a later time, one handedness 
is predominant in the system. 

1.1 Models for homochiralisation 

Many models have been proposed for the emergence of homochirality from an initially racemic 
mixture of precursors. Frank [10] proposed an open system into which R and S particles are 
continually introduced, and combine to form one of two possible products: left- or right-handed 



species, X,Y. Each of these products acts as a catalyst for its own production (autocatalysis), 
and each combines with the opposing handed product (cross-inhibition) to form an inert product 
(P) which is removed from the system at some rate. These processes are summarised by the 
following reaction scheme: 

external source — Y R, S input, k , 

R + S ^ X R + S # Y slow, ^, 

R + S + X t± 2X R + S + Y # 2Y fast, auto catalytic, k 2 (1.1) 

X + y — > P cross-inhibition, fc 3 , 

P — y removal, k 4 . 

Ignoring the reversible reactions (for simplicity), this system can be modelled by the differential 
equations 

— = k - 2kirs - k 2 rs(x+y) + k-i(x+y) + k^ 2 (x 2 +y 2 ), (1.2) 

— = k - 2k]rs- k 2 rs(x+y) + k-i(x+y) + k_ 2 (x 2 +y 2 ), (1.3) 

— — = k\rs + k 2 rsx — k^xy — k_\x — k^ 2 x 2 , (1-4) 

at 

— = k ± rs + k 2 rsy - k 3 xy - k^ ± y - k„ 2 y 2 , (1.5) 

— = k 3 xy-k 4 p, (1.6) 

from which we note that at steady-state we have 

k + k_ t (x + y) + k^(x 2 + y 2 ) 

rs = ; T-, n • (1-7 ) 

2^ + hix + y) V ; 

We write the absolute enantiomeric excess as ee = x — y and the total concentration as a = x + y; 
adding and subtracting the equations for dx/dt and dy/dt, we find 

a = bee , (1.8) 
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k 2 (k_ 2 ee 2 + k_ 2 a 2 + 2k_ l( T + 2k ) , 

k-i — k_ 2 a 



0. (1.9) 



2(2*4 + k 2 a) 

Hence ee = is always a solution, and there are other solutions with ee ^ if the rate constants 
k* satisfy certain conditions (these include k% > k_ 2 and ko being sufficiently large). 
The important issues to note here are: 

(i) this system is open, it requires the continual supply of fresh R, S to maintain the asymmetric 
steady-state. Also, the removal of products is required to avoid the input terms causing 
the total amount of material to increase indefinitely; 

(ii) the forcing input term drives the system away from an equilibrium solution, into a distinct 
steady-state solution; 

(iii) the system has cross-inhibition which removes equal numbers of X and Y, amplifying any 
differences caused by random fluctuations in the initial data or in the input rates. 



Saito & Hyuga [23] discuss a sequence of toy models describing homochirality caused by 
nonlinear autocatalysis and recycling. Their family of models can be summarised by 

— = kr 2 (l - r - s) - Ar, (1.10) 

ds 

=- = ks 2 {l-r-s)-Xs, (1.11) 

where r and s are the concentrations of the two enantiomers. Initially they consider k r = k s = k 
and A = and find that enantiomeric exess, r — s is constant. Next the case k r = kr, k s = ks, 
A = is analysed, wherein the relative enantiomeric excess ^rf is constant. Then the more 
complex case of k r = kr 2 , k s = ks 2 , A = is analysed, and amplification of the enantiomeric 
excess is obtained. This amplification persists when the case A > is finally analysed. This 
shows us strong autocatalysis may cause homochiralisation, but in any given experiment, it is 
not clear which form of rate coefficients (k r , k s , A) should be used. 

Saito & Hyuga (2005) analyse a series of models of crystallisation which include some of 
features present in our more general model. They note that a model truncated at tetramers 
exhibits different behaviour from one truncated at hexamers. In particular, the symmetry- 
breaking phenomena is not present in the tetramer model, but is exhibited by the hexamer 
model. Hence, later, we will consider models truncated at the tetramer and the hexamer levels 
and investigate the differences in symmetry-breaking behaviour (Sections |3] and H|). 

Denoting monomers by c, small and large left-handed clusters by X\,x 2 respectively and 
right-handed by 2/1,2/2, Uwaha [2E] writes down the scheme 

dc 

— = -2k z 2 k 1 z(x 1 + 2/1) + Xi(x 2 + y 2 ) + A (xi + yi), (1.12) 

at 

— — = k z 2 - k u x 1 x 2 - k c x\ + X u x 2 + A xi, (1-13) 

at 

— kix 2 c + k u x\x 2 + k c x\ - \\x 2 - X u x 2 , (1-14) 



dt 

k z 2 - k u y ± y 2 - k c y 2 + \ u y 2 + \ y 1 , (1.15) 



d ?/l _ ,. 2 , , -2 



dt 

-— = k x y 2 c + k u yxy 2 + k c y\ - X 1 y 2 - A u j/ 2 , (1.16) 

at 

which models 

• the formation of small chiral clusters (xi,yi) from an achiral monomer (c) at rate ko, 

• small chiral clusters (xi,yi) of the same handedness combining to form larger chiral clusters 
(rate k c ), 

• small and larger clusters combining to form larger clusters (rate k u ), 

• large clusters combining with achiral monomers to form more large clusters at the rate ki, 

• the break up of larger clusters into smaller clusters (rate X u ), 

• the break up of small clusters into achiral monomers (rate A ), 

• the break up of larger clusters into achiral monomers (rate Ai). 



Such a model can exhibit symmetry-breaking to a solution in which x\ 7^ X2 and X2 7^ yi- 
Uwaha points out that the recycling part of the model (the A* parameters) are crucial to the 
formation of a 'completely' homochiral state. One problem with such a model is that since the 
variables are all total masses in the system, the size of clusters is not explicitly included. In 
asymmetric distributions, the typical size of left- and right- handed clusters may differ drastically, 
hence the rates of reactions will proceed differently in the cases of a few large crystals or many 
smaller crystals. 

Sandars has proposed a model of symmetry-breaking in the formation of chiral polymers 
[25] . His model has an achiral substrate (S) which splits into chiral monomers Li,Ri both 
spontaneously at a slow rate and at a faster rate, when catalysed by the presence of long 
homochiral chains. This catalytic effect has both autocatalytic and crosscatalytic components, 
that is, for example, the presence of long right-handed chains R n autocatalyses the production 
of right-handed monomers Ri from S, (autocatalysis) as well as the production of left-handed 
monomers, L\ (crosscatalysis). Sandars assumes the growth rates of chains are linear and not 
catalysed; the other mechanism required to produce a symmetry-breaking bifurcation to a chiral 
state is cross-inhibition, by which chains of opposite handednesses interact and prevent either 
from further growth. These mechanisms are summarised by 

S — > Li, S — V Ri, slow, 

S+L n — > L\ + L n , S+R n ^ Ri + R n , autocatalytic, rate oc 1 + /, 

S+R n — > Li + R n , S+L n — y Ri+L n , cross-catalytic, rate oc 1 — /, 

L n + L\ — \ L n+ i, R n + R\ — > R n+ i, chain growth, rate = a, 

L n + R\ — > Qn+i, R n + Li — > P n +i, cross- inhibition, rate = a\- 

This model and generalisations of it have been analysed by Sandars [25], Brandenburg et al. [HIE], 
Multimaki & Brandenburg [19], Wattis & Coveney [32], [33], Gleiser & Walker [11], Gleiser et al. 
[T2] . Typically a classic pitchfork bifurcation is found when the fidelity (/) of the autocatalysis 
over the cross-catalysis is increased. One counterintuitive effect is that increasing the cross- 
inhibition effect (x) aids the bifurcation, allowing it to occur at lower values of the fidelity 
parameter /. 

1.2 Experimental results on homo chiralisat ion 

The Soai reaction was one of the first experiments which demonstrated that a chemical reaction 
could amplify initial small imbalances in chiral balance; that is, a small enantiomeric exess in 
catalyst at the start of the experiment led to a much larger imbalance in the chiralities of the 
products at the end of the reaction. Soai et al. [27] was able to achieve an enantiomeric exess 
exceeding 85% in the asymmetric autocatalysis of chiral pyrimidyl alkanol. 

The first work showing that crystallisation experiments could exhibit symmetry breaking was 
that of Kondepudi & Nelson [T5] . Later Kondepudi et al. [2] showed that the stirring rate was 
a good bifurcation parameter to analyse the final distribution of chiralities of crystals emerging 
from a supersaturated solution of sodium chlorate. With no stirring, there were approximately 
equal numbers of left- and right-handed crystals. Above a critical (threshold) stirring rate, the 
imbalance in the numbers of each handedness increased, until, at large enough stirring rates, total 
chiral purity was achieved. This is due to all crystals in the system being derived from the same 
'mother' crystal, which is the first crystal to become established in the system; all other crystals 
grow from fragments removed from it (either directly or indirectly). Before this, Kondepudi 
& Nelson [T6| [T7] worked on the theory of chiral symmetry-breaking mechanisms with the aim 



of predicting how parity-violating perturbations could be amplified to give an enantiomeric 
exess in prebiotic chemistry, and the timescales involved. Their results suggest a timescale of 
approximately 10 4 years. More recently, Kondepudi and Asakura [T3] have summarised both the 
experimental and theoretical aspects of this work. 

Viedma |29j was the first to observe that grinding a mixture of chiral crystals eventually led 
to a distribution of crystals which were all of the same handedness. The crystalline material 
used was sodium chlorate, as used by Kondepudi et al. [15J. Samples of L and D crystals are 
mixed with water in round-bottomed flasks and the system is stirred by a magnetic bar (of 
length 3-20mm) at 600rpm. The system is maintained in a supersaturated state; small glass 
balls are added to continually crush the crystals. The grinding is thus continuous, and crystals 
are maintained below a size of 200 /im. The chirality of the resulting crystals was determined by 
removing them from the flask, allowing them to grow and measuring their optical activity. The 
results show that, over time, the percentages of left- and right-handed crystals steadily change 
from about 50/50 to 100/0 or 0/100 - a state which is described as complete chiral purity. With 
stirring only and no glass balls, the systems conserve their initial chiral excesses; with glass balls 
present and stirring, the chiral excess increases, and this occurs more rapidly if more balls are 
present or the speed of stirring is increased. 

More recently, Noorduin et al. [2T] have observed a similar effect with amino acids - a much 
more relevant molecule in the study of origins of life. This work has been reviewed by McBride 
& Tully [18], who add to the speculation on the mechanisms responsible for the phenomenon. 
Noorduin et al. describe grinding as 'dynamic dissolution/crystallization processes that result 
in the conversion of one solid enantiomorph into the other'. They also note that 'once a state 
of single chirality is achieved, the system is "locked" because primary nucleation to form and 
sustain new crystals from the opposite enantiomer is kinetically prohibited'. Both these quotes 
include the crucial fact that the process evolves not towards an equilibrium solution (which 
would be racemic), but towards a different, dynamic steady-state solution. As noted by Plasson 
(personal communication, 2008), this nonequilibrium state is maintained due to the constant 
input of energy into the system through the grinding process. 

McBride & Tully [IB] discuss the growth of one enantiomorph, and the dissolution of the other 
as a type of Ostwald ripening process; with the large surface area to volume ratio of smaller 
crystals giving a rapid dissolution rate, whilst larger crystals, have a lower surface area to volume 
ratio meaning that they dissolve more slowly. However appealing such an argument maybe, since 
surface area arguments can equally well be applied to the growth side of the process, it is not 
clear that this is either necessary or sufficient. Infact, the model analysed later in this paper will 
show that a critical cluster size is not necessary to explain homochiralisation through grinding. 

1.3 Our aims 

We aim to describe the results of the crystal grinding phenomenon through a model which recy- 
cles mass through grinding, which causes crystals to fragment, rather than having explicit mass 
input and removal. Simultaneously we need crystal growth processes to maintain a distribution 
of sizeable crystals. 

We assume that the crystals are solids formed in an aqueous environment, however, we leave 
open questions as to whether they are crystals of some mineral of direct biological relevance 
(such as amino acids), or whether they are some other material, which after growing, will later 
provide a chirally selective surface for biomolecules to crystallise on, or be a catalyst for chiral 
polymerisation to occur. Following Darwin's [9] "warm little pond" , an attractive scenario might 
be a tidal rock pool, where waves agitating pebbles provide the energetic input for grinding. 



Taking more account of recent work, a more likely place is a suboceanic hydrothermal vent 
where the rapid convection of hot water impels growing nucleii into the vent's rough walls as 
well as breaking particles off the walls and entraining them into the fluid flow, simultaneously 
grinding any growing crystals. 

In Section [2] we propose a detailed microscopic model of the nucleation and crystal growth 
of several species simultaneously This has the form of a generalised Becker-Doring system of 
equations pp. Due to the complexity of the model we immediately simplify it, making assump- 
tions on the rate coefficients. Furthermore, to elucidate those processes which are responsible 
for homo chiralisat ion, we remove some processes completely so as to obtain a simple system of 
ordinary differential equations which can be analysed theoretically. 

The simplest model which might be expected to show homochiralisation is one which has 
small and large clusters of each handedness. Such a truncated model is considered in Section 
[3] wherein it is shown that such a model might lead to amplification of enantiomeric exess in 
the short time, but that in the long-time limit, only the racemic state can be approached. This 
model has the structure akin to that of Saito & Hyuga [21] truncated at the tetramer level. 

Hence, in Section H] we consider a more complex model with a cut-off at larger sizes (one 
can think of small, medium, and large clusters of each handedness). Such a model has a similar 
structure to the hexamer truncation analysed by Saito & Hyuga [21] . We find that such a model 
does allow a final steady-state in which one chirality dominates the system and the other is 
present only in vanishingly small amounts. 

However, as discussed earlier, there may be subtle effects whereby it is not just the number 
of crystals of each type that is important to the effect, but a combination of size and number of 
each handedness of crystal that is important to the evolution of the process. Hence, in Section 
we introduce an alternative reduction of the system of governing equations. In this, instead 
of truncating and keeping only clusters of a small size, we postulate a form for the distribution 
which includes information on both the number and size of crystals, and use these two quantities 
to construct a system of five ordinary differential equations for the system's evolution. 

We discuss the results in Sections |6] and [7] which conclude the paper. The Appendix IA1 shows 
how, by removing the symmetry in the growth rates of the two handednesses, the model could 
be generalised to account for the competitive nucleation of different polymorphs growing from a 
common supply of monomer. 

2 The BD model with dimer interactions and an amor- 
phous metastable phase 

2.1 Preliminaries 



Smoluchowski [2S] proposed a model in which clusters of any sizes could combine pairwise to 
form larger clusters. Chemically this process is written C r + C s — )■ C r+S where C r represents a 
cluster of size r. Assuming this process is reversible and occurs with a forward rate given by a TyS 
and a reverse rate given by b r>s , the law of mass action yields the kinetic equations 



dc r r ~ 1 °° 

— — — — — y ^ \Qj sr — s C s C r — s sr — s C r ) y j \U>r,sCrCs rs C r 
at s=l s=l 



(2.1) 



These are known as the coagulation-fragmentation equations. There are simplifications in which 
only interactions between clusters of particular sizes are permitted to occur, for example when 

6 



only cluster-monomer interactions can occur, the Becker-Doring equations pQ are obtained. Da 
Costa has formulated a system in which only clusters upto a certain size (N) are permitted to 
coalesce with or fragment from other clusters. In the case of N = 2, which is pertinent to the 
current study, only cluster-monomer and cluster-dimer interactions are allowed, for example 

C r + Ci ■*— G r +i, C r + C2 -r- w+2- (2-2) 

This leads to a system of kinetic equations of the form 

J r _ ± - J r + K r _ 2 - K r , (r>3), (2.3) 

00 
J1-J2- tf 2 -£#,., (2.4) 



dc r 

~dt 

dc 2 
~dt 
dci 

"dT 



00 



Ji-K 2 -Y,Jr, (2.5) 



r=l 



J r — a r C r C\ — 6 r _)_iC r _)_l, i\ r — a r C r C2 — Pr+2 c -r+2- (2-6) 

A simple example of such a system has been analysed previously by Bolton & Wattis [3] . 

In the next subsection we generalise the model (12.11) to include a variety of 'species' or 'mor- 
phologies' of cluster, representing left-handed, right-handed and achiral clusters. We simplify 
the model in stages to one in which only monomer and dimer interactions are described, and 
then one in which only dimer interactions occur. 

2.2 A full microscopic model of chiral crystallisation 

We start by outlining all the possible cluster growth, fragmentation and transformation pro- 
cesses. We denote the two handed clusters by X r , Y r , where the subscript r specifies the size of 
cluster. Achiral clusters are denoted by C r , and we allow clusters to change their morphology 
spontaneously according to 

(2-7) 

We allow clusters to grow by coalescing with clusters of similar handedness or an achiral cluster. 
In the case of the latter process, we assume that the cluster produced is chiral with the same 
chirality as the parent. Thus 



(2-8) 



We do not permit clusters of opposite to chirality to merge. Finally we describe fragmentation: 

all clusters may fragment, producing two smaller clusters each of the same chirality as the parent 

cluster 

X r + S — y X r -\- X s rate = p r s , 

o^-j-s t Oy ~t~ o s rate — — 6^ s , (z.yj 

±r-\-s ' *r ' * s rate — — p r s . 



, - 


-». x r 


TcltG — l^V) 


Jv r 


-)• c r 


Fclte — yb r v r ^ 


C r 


-»- Y r 


TcltG == l^V) 


Y r 


-)- C r 


rate = fi r u r . 



X r + X s 


-)■ 


X r + S , 


rate = 


= sr,s; 


X r -\- G s 


-)■ 




rate = 


~ "r,S) 


O r T" O s 


— 7 


c 

KJ r+si 


rate = 


~ ^r,si 


Y r + C s 


— » 


v 

± r+S) 


rate = 


~ ®-r,si 


Y r + Y s 


->■ 


1 r+si 


rate = 


~ Kr,s- 



Setting up concentration variables for each size and each type of cluster by defining c r (t) = [C r ], 
x r (t) = [X r ], y r (t) = [Y r ] and applying the law of mass action, we obtain 



dt 



dt 



dc r °° 

— = -2n r Cr + /j, r v r (x r + y r ) - ^2 a kt rC r (x k + Vk) (2.10) 

fc=l 

1 — 1 oo 

+\ E (h,r- k c k c r . k - e k , r . k c k c r . k ) - £ (5 k , r c k c r - e Kr c r+k ) , 
fe=i fc=i 

r— 1 r— 1 

H r C r fi r V r X r t / Oi k ^ r — k C k Xj — fc 2 / , (.Sfc,? - — k%k%l — fc Pk,r—k%r ) 
k=l fe=l 

oo 

— 2_^ \Sk,r x k x r ~ Pk,r x r+k) j (2-H) 

fc=l 

dy r r_1 r_1 

— — = ll r C T — [l r V r y r + E ®k,r-kCkyr-k—2 E (6c,r-fc2/fc2/r-fc ~ Pk,r-kVr) 
at fc=l fc=l 

oo 

- E (£k, r yky r - (3k, r y r +k) • (2.12) 

fc=i 

The main problem with such a model is the vast number of parameters that have been introduced 

(ttr.fc, £r,fc) /?r,fe; /^r, ^V> ^r,fc> e r,k, for all fc,r). 

Hence we make several simplifications: 

(i) we assume that the dominant coagulation and fragmentation processes are between large 
and very small clusters (rather than large clusters and other large clusters). Specifically, 
we assume that only coalescences involving C\ and Ci need to be retained in the model, 
and fragmentation always yields either a monomer or a dimer fragment. This assumption 
means that the system can be reduced to a generalised Becker-Doring equation closer to 
the form of Q-Q rather than (TO) : 

(ii) we also assume that the achiral clusters are unstable at larger size, so that their presence is 
only relevant at small sizes. Typically at small sizes, clusters are amorphous and do not 
take on the properties of the bulk phase, hence at small sizes clusters can be considered 
achiral. We assume that there is a regime of cluster sizes where there is a transition to 
chiral structures, and where clusters can take on the bulk structure (which is chiral) as 
well as exist in amorphous form. At even larger sizes, we assume that only the chiral forms 
exist, and no achiral structure can be adopted; 

(iv) furthermore, we assume that all rates are independent of cluster size, specifically, 

(r > 2) (2.13) 

(2.14) 
(2.15) 
(2.16) 
(2.17) 
(2.18) 



a k,l 


= a, 




U fc2 ^' ^k,r "5 


/'2 


= V, 




l_i r = 0, (r > 3), 


y-i 


= ", 




v r = 0, (r > 3), 


*L,1 


= 5, 


Ok,r 


= 0, (otherwise) 


e l,l 


= e, 


€-k,r 


= 0, (otherwise) 


£fc,2 


— &,k : 


= £, 


6s,r = 0, (otherwise) 


0k,l 


= Pl,k 


= b, 


Pk,2 = @2,k = P, Pk,r 



0, (otherwise), 



(2.19) 



Ultimately we will set a = b = = 5 = e so that we have only five parameters to consider 
(a, f, 0, fi, v). 
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Figure 1: Reaction scheme involving monomer and dimer aggregation and fragmentation of 
achiral clusters and those of both handednesses (right and left). The aggregation of achiral and 
chiral clusters is not shown (rates a, £). 



This scheme is illustrated in Figure [TJ However, before writing down a further system of 
equations, we make one further simplification. We take the transition region described in (ii), 
above, to be just the dimers. Thus the only types of achiral cluster are the monomer and the 
dimer (ci, c 2 ); dimers exist in achiral, right- and left-handed forms (c 2 , X2, 1/2); at larger sizes 
only left- and right-handed clusters exist (x r , y r , r > 2). 

The kinetic equations can be reduced to 



dci 
"dt" 

dc 2 
~dt 
dx r 
~dT 

dx 3 

dt 
dx 2 
~df 



dy r 
dt 

dj/3 
dt 

dt 



2ec2 — 25c\ — ^2(acix r + ac\y T — bx r+ \ — by r+ i), 

r=2 

00 

8c\ - SC 2 - 2/iC 2 + fJLv(x2 + 2/2) - ^ aC z( X r + Vr), 

r=2 

ac\x r -\ — bx r — ac\x r + bx r+ \ + ac2X r -2 — ac2X r 
-ftx r + ftx r+2 + ix 2 x r ^2 - ix 2 x ri (r > 4), 
ac\X2 — bx 3 — CLC1X3 + 6x4 — ac 2 x 3 — C,x 2 x 3 + ftx 5 , 

/iC2 — 11VX2 + bx% — ac\X2 — CXX2C2 + ftxi 
00 00 

+ Yl ft X r+2 ~ J2 &2 X r ~ &2, 

r=2 r=2 

ac x y r -\ - by r - ac x y r + by r+1 + ac 2 y r ~2 - ac 2 y r 

-ftVr + ftyr+2 + &/22/r-2 ~ ^WlVr, (?" > 4), 

ocij/2 - by 3 - aciy 3 + by 4 - ac 2 y 3 - ^y 2 ys + ftVa, 
fxc 2 - nvy 2 + by 3 - ac t y 2 - ay 2 c 2 + fty 4 



+ J2 ftvr+2 - J2 tv*ur - €v 



(2.20) 
(2.21) 

(2.22) 
(2.23) 

(2.24) 

(2.25) 
(2.26) 

(2.27) 



r=2 



r=2 



2.3 Summary and simulations of the macroscopic model 

The advantage of the above simplifications is that certain sums appear repeatedly; by defining 
new quantities as these sums, the system can be written in a simpler fashion. We define N x = 
2 x r , N y = Y,T=2 Vr, then 



/ v i 9 ^Vi -t'fl 



U.09 o 

-j— = (^Cj -ec 2 - 2//c 2 + ^i/(x2 + 2/2) — cuc 2 (N x + iV 2 



dci 

It" 

dc 2 

It 

dt 
dx 2 

dNy_ 

dt 

dy 2 

dt 



2ec 2 - 25c? - aci(iV x + N y ) + b{N x - x 2 + N y - y 2 ) 



yj> 



fic 2 - \ivx 2 + 0(N X - x 3 - x 2 ) - £x 2 N x , 

\xc 2 — nux 2 + bx 3 — ac\X 2 — ax 2 c 2 + fi(x± + N x — x 2 — x 3 ) 

-£x 2 2 - £x 2 N x , 

Hc 2 - [ivy 2 + /3(N y -y 3 - y 2 ) - £y 2 N y , 

\ic 2 - ixvy 2 + by 3 - ac x y 2 - ay 2 c 2 + f3(y 4 + N y - y 2 - y 3 ) 



(2.28) 
(2.29) 
(2.30) 

(2.31) 
(2.32) 

(2.33) 



However, such a system of equations is not 'closed'. The equations contain x 3 , y 3 , x 4 , j/ 4 , and yet 
we have no expressions for these; reintroducing equations for x 3 ,y 3 would introduce £5,3/5 and 
so an infinite regression would be entered into. 
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Figure 2: Plot of the concentrations c%, c 2 , N x , N y , N = N x + N y , g x , g y , g x + g y and 
Qx + By + 2c 2 + cl against time, t on a logarithmic timescale. Since model equations are in 
nondimensional form, the time units are arbitrary. Parameter values fi = 1.0, v = 0.5, 5 = 1, 
e = 5, a = 4, b = 0.02, a = 10, f = 10, fi = 0.03, with initial conditions c 2 = 0.49, x 4 (0) = 0.004, 
2/4(0) = 0.006, and all other concentrations zero. 

Hence we need to find some suitable alternative expressions for x 3 , y 3 , £4, y 4 ; or an alternative 
way of reducing the system to just a few ordinary differential equations that can easily be 
analysed. Such systems are considered in Sections El IH and El Before that, however, we illustrate 
the behaviour of the system by briefly presenting the results of some numerical simulations. In 
Figures [2] and [3] we show the results of a simulation of f)2.28l) - fl2.33l) . The former shows the 
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Figure 3: Plot of the cluster size distribution at t = (dashed line), t — 112 (dotted line) and 
t = 9.4 x 10 5 . Parameters and initial conditions as in Figure EJ 

evolution of the concentrations c\ which rises then decays, Ci which decays since the parameters 
have been chosen to reflect a cluster-dominated system. Also plotted are the numbers of clusters 
N x , N y and the mass of material in clusters g x , g y defined by 



A' 



A 



i=2 



X 



3i 



Qy 



Y^JVi- 



(2.34) 



i=2 



Note that under this definition g x + g y + c\ + 2c 2 is conserved, and this is plotted as rho. Both 
the total number of clusters, N x + N y , and total mass of material in handed clusters q x + g y 
appear to equilibrate by t — 10 2 , however, at a much later time (t ~ 10 4 — 10 5 ) a symmetry- 
breaking bifurcation occurs, and the system changes from almost racemic (that is, symmetric) 
to asymmetric. This is more clearly seen in Figure El where we plot the cluster size distribution 
at three time points. At t = there are only dimers present (dashed line), and we impose a 
small difference in the concentrations of X2 and y 2 - At a later time, t = 112 (dotted line), there is 
almost no difference between the X- and ^-distributions, however by the end of the simulation 
(t ~ 10 6 , solid line) one distribution clearly completely dominates the other. 



2.4 Simplified macroscopic model 

To obtain the simplest model which involves three polymorphs corresponding to right-handed and 
left-handed chiral clusters and achiral clusters, we now aim to simplify the processes of cluster 
aggregation and fragmentation in (I2.28l) - fl2.33l) . Our aim is to retain the symmetry-breaking 
phenomenon but eliminate physical processes which are not necessary for it to occur. 

Our first simplification is to remove all clusters of odd size from the model, and just consider 
dimers, tetramers, hexamers, etc. This corresponds to putting a = 0, b = which removes 
X3 and 1/3 from the system. Furthermore, we put e = and make S large, so that the achiral 
monomer is rapidly and irreversibly converted to achiral dimer. Since the monomers do not then 
influence the evolution of any of the other variables, we further simplify the system by ignoring 
C\ (or, more simply, just impose initial data in which Ci(0) = 0). Thus we are left with 



dc 2 

"dT 



-2fj,c 2 + \iv{x 2 + 1/2) - ctc 2 {N x + N t 



V)i 



(2.35) 
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dN x 

dt 
dx 2 

~dT 

dNy 

dt 

dt 



fic 2 - \xvx 2 + /3{N X - x 2 ) - £x 2 N x , 

fic 2 - [ivx 2 - ax 2 c 2 + /3(N x -x 2 + x A ) - £,xj - C,x 2 N x , 

/j,c 2 - fiuy 2 + /3(N y - y 2 ) - £y 2 N y , 

Hc 2 - jj,uy 2 - ay 2 c 2 + /3(N y -y 2 +y 4 ) - ^y\ - £y 2 N y . 



(2.36) 
(2.37) 
(2.38) 
(2.39) 



Since we have removed four parameters from the model, and halved the number of dependent 
variables, we show a couple of numerical simulations just to show that the system above does 
still exhibit symmetry-breaking behaviour. 
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Figure 4: Plot of the concentrations c%, c 2 , N x , N y , N 



N* 



N y , Q x , Qy, Q x + Qy and 
Qx + By + 2c2 + C\ against time, t on a logarithmic timescale. Since model equations are in 
nondimensional form, the time units are arbitrary. Parameter values fi = I, v = 0.5, a = 10, 
£ = 10, (3 = 0.03, with initial conditions c 2 = 0.49, x 4 (0) = 0.004, y 4 (0) = 0.006, all other 
concentrations zero. 

Figure H] appears similar to Figure |2l suggesting that removing the monomer interactions 
has changed the underlying dynamics little. We still observe the characteristic equilibration of 
cluster numbers and cluster masses as c 2 decays, and then a period of quiesence (t ~ 10 to 10 4 ) 
before a later symmetry-breaking event, around t ~ 10 5 . At first sight, the distribution of X- 
and ^-clusters displayed in Figure [5] is quite different to Figure [3j this is due to the absence of 
monomers from the system, meaning that only even-sized clusters can now be formed. If one 
only looks at the even-sized clusters in Figure O we once again see only a slight difference at 
t = (dashed line), almost no difference at t « 250 (dotted line) but a significant difference at 
t = 6 x 10 5 (solid line). We include one further graph here, Figure M similar to Figure H] but 
on a linear rather than a logarithmic timescale. This should be compared with Figures such as 
Figures 3 and 4 of Viedma [29J and Figure 1 of Noorduin et al. [21]. 
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Figure 5: Plot of the cluster size distribution at t = (dashed line), t = 250 (dotted line) and 
t = 6 x 10 5 . Parameters and initial conditions as in Figure HI 
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Figure 6: Plot of the concentrations C\, C2, N x , N y , N = N x + N y , g x , g y , g x + g y and 
Qx + Qy + 2c2 + c\ against time, t on a logarithmic timescale. Parameters and initial conditions 
as in Figure HI 
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Figure 7: Simplest possible reaction scheme which might exhibit chiral symmetry-breaking. 
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3 The truncation at tetramers 

The simplest possible reaction scheme of the form ( 12.201) — ( 0.27J) which we might expect to exhibit 
symmetry-breaking to homochirality is the system truncated at tetramers, namely 

dco 

— = -2/jc 2 + fiis(x 2 + y 2 ) -ac 2 (x 2 +V2J, (3.1) 

— /ic 2 — \ivx 2 — ac 2 x 2 — 2^x\ + 2/3x4, (3.2) 



dt 

/j,c 2 - fii/y 2 - ac 2 y 2 - 2^y\ + 2/% 4 , (3.3) 



d ?/2 _ ,,,.;! 



dt 

—- ax 2 c 2 + ^x\ — (3x4, (3.4) 



df 

ay 2 c 2 + £y 2 2 - f3y 4 . (3.5) 



d?/ 4 , 2 



to 



dt 

We investigate the symmetry-breaking by transforming the variables x 2 , x 4 , y 2 , 7/4 according 

x 2 = \z(l + 9), y 2 =\z{\-9), (3.6) 

X4,= \w(l + 4>), y A = \w(l- 4>), (3.7) 

where z = x 2 + y 2 is the total concentration of chiral dimers, w = x<± + 1/4 is the total tetramer 
concentration, 9 = (x 2 — y 2 )/z is the relative chirality of the dimers, = (24 — y±)/w is the 
relative chirality of tetramers. Hence 

— — = — 2fic 2 + fiuz — ac 2 z, (3.8) 

C.1 v 

dz 

— = 2fic 2 - ^uz-ac 2 z-^z 2 (l + 9 2 ) + 2(3w, (3.9) 

dt 

^ = azc 2 + ^z 2 (l + 9 2 )-Pw, (3.10) 

d i = 4fc + ^ + ^-^^, (3.11) 

dt \ z z ) z 

^ = e-(ac + Zz)-(ac+Uz(l + 9 2 ))-<p. (3.12) 

dt w v 2 ' w 

The stability of the evolving symmetric-state (9 = <fi = 0) is given by the eigenvalues (q) of the 
matrix 

'■241c _|_ 2§w_ _|_ t^ M^ 

(ac + ^)^ -(ac+±£z) 



(3.13) 



which are given by 



9 1 acz £z 2/ic . 2/3um 
q 2 + q[ + — + — + &+ -^— + 

\ W W Z Z J 



— (2fj,cac + [ic£z + ac^z 2 + ±£V - p£zw) = 0. (3.14) 

w \ * ) 

Hence there is an instability if 

/3£zw > 2ficac + fic^z + ac£z 2 + f £ 2 £ 3 , (3.15) 
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using the steady-state result that 2(3w = z(2ac + £z) and factorising (2ac + £z) out of the 
result, reduces the instability (13. 15[) to the contradictory ^z 2 > ^z 2 + 2/ic. Hence the racemic 
steady-state of the system is stable for all choices of parameter values and is approached from 
all initial conditions. However, initial perturbations, may be amplified due to the presence of 
nonlinear terms. 
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Figure 8: The concentrations C2, z and w f l3.6p -( |3T7l) plotted against time, for the tetramer- 
truncated system with the two sets of initial data (I3.16p . Since model equations are in nondimen- 
sional form, the time units are arbitrary. The parameter values are // = 1, v = 0.5, a = £ = 10, 
/3 = 0.1. 
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Figure 9: The chiralities 6, <p (l3.6p -f l3~7T|) plotted against time, for the tetramer-truncated system 
with the two sets of initial data (I3.16p . Since model equations are in nondimensional form, the 
time units are arbitrary. The parameter values are the same as in Figure [8j 



Evolution from two sets of initial conditions of the system (I3.ip - fl3.5p are shown in each of 
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Figures [HI |9j The continuous and dotted lines correspond to the initial data 

c 2 (0) = 0.29, x 2 (0) = 0.0051, y 2 {0) = 0.0049, 

x 4 (0) = 0.051, ty 4 (0) = 0.049; and . , 

c 2 (0) = 0, x 2 (0) = 0.051 y 2 (0) = 0.049, [ ' 

x 4 (0) = 0.1, y 4 (0) = 0.1; 

respectively. In the former case, the system starts with considerable amount of amorphous dimer, 
which is converted into clusters, and initially there is a slight chiral imbalance in favour of x 2 
and £ 4 over y 2 and 7/4. Over time this imbalance reduces (see figure |H]); although there is a region 
around t = 1 where 9 increases, both 9 and eventually approach the zero steady-state. 

For both sets of initial conditions we note that the chiralities evolve over a significantly longer 
timescale than the concentrations, the latter having reached steady-state before t = 10 and the 
former still evolving when t = O(10 2 ). In the second set of initial data, there is no c 2 present 
initially and there are exactly equal numbers of the two chiral forms of the larger cluster, but a 
slight exess of x 2 over y 2 . In time an imbalance in larger clusters is produced, but over larger 
timescales, both 9 and <p again approach the zero steady-state. 

Hence, we observe that the truncated system ( I3.ip - (j3.5p does not yield a chirally asymmet- 
ric steady-state. Even though in the early stages of the reaction chiral perturbations may be 
amplified, at the end of the reaction there is a slower timescale over which the system returns to 
a racemic state. In the next section we consider a system truncated at hexamers to investigate 
whether that system allows symmetry-breaking of the steady-state. 

4 The truncation at hexamers 

The above analysis has shown that the truncation of the model ( I2.20p -( ]2.27p to ( I3.ip -( l3~5|) 
results in a model which always ultimately approaches the symmetric (racemic) steady-state. In 
this section, we show that a more complex model, the truncation at hexamers retains enough 
complexity to demonstrate the symmetry-breaking bifurcation which occurs in the full system. 
In this case the governing equations are 

dc 2 

"d7 

dx 2 _ 0f ., 

~df 

dx4 2 



dx 6 

~df 

dj/2 

dt 



-2/ic 2 + /!z/(x 2 + 2/2) - ac 2 (x 2 + 2/2) - ac 2 (xi + 7/4), (4.1) 

[ic 2 — [ivx 2 — ac 2 x 2 — 2£x 2 — £x 2 X4 + 2/3x4 + (3xq, (4.2) 

ax 2 c 2 + £x 2 — /3x4 — ac 2 X4 — £,x 2 Xi + (3xq, (4.3) 

ax4C 2 + £x 2 x 4 — [3x6, (4-4) 

jj,c 2 - fj,i/y 2 - ac 2 y 2 - 2£yl - ^y 2 y 4 + 2/3y 4 + f3y 6 , (4.5) 



— = ay 2 c 2 + iy\ - (3y 4 - ac 2 y 4 - £y 2 y 4 + /3y e , (4.6) 

_ 1 = ay 4 c 2 + iy 2 y 4 - /3y 6 . (4.7) 

To analyse the symmetry-breaking in the system we transform the dependent coordinates 
from x 2 ,X4,XQ,y 2 ,y4,y 6 to total concentrations z,w,u and relative chiralities 9,4>,ip according 

to 

x 2 = ~z(l + 9), x A = |ty(H-0), xq = ~u(l + ip), 

(4.8) 
y 2 = \z{l-9), y 4 = |ty(l-0), y 6 = |u(l -?/>). 
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We now separate the governing equations for the total concentrations of dimers (c, z), tetramers 
(w) and hexamers (u) 

dc 

— = — 2pLC + fivz — acz — acw, (4.9) 

CI / 

dz 

— = 2/ic- fivz-acz-£z 2 (l + 9 2 ) - \zw(\ + 0<f>) + (3u + 2f3w, 

kJLL 



dw 

~dt 
du 

dt 

from those for the chiralities 



(4.10) 
acz + \£z 2 (l + d 2 ) -flw + pu- acw - \&w(\ + 9(f)), (4.11) 

acw + \£zw(l + 6(f)) - (3u, (4.12) 



£ = azu-Q + epLv+t-t-m ( 4.i3) 

dt u lu 

^ = — (e-4>) + ^(29 -0-^) + -(V> _^)_i^(i-A 

dt w 2w w 



d6 2/xcfl 2 y 2 (3wp f3u6 

2/3w0 2/3w# 



(4.14) 



(4.15) 



z z 

In applications, we expect v < 1, so that the small amorphous clusters (dimers) prefer to 
adopt one of their chiral states rather than the achiral structure. In addition, we note that the 
grinding process observed in experiments is much longer than the crystallisation process, and 
that there are many larger, macroscopic crystals hence we consider two limits in which j3 <C a£. 
We will consider the case of small /3 with all other parameters being 0(1) and then the case 
where a ~ £ 3> 1 and all other parameters are 0(1). 

4.1 Symmetric steady-state for the concentrations 

Firstly, let us solve for the symmetric steady-state. In this case we assume 9 = = <fi = ip, 
simplifying equations ( l4.9p -( l4"T2"j) . One of these is a redundant equation, hence we have the 

solution 

z z 

w = -(ac+ z&), u= —(ac+ \£z) 2 , (4.16) 

f.|f£ + ftf + ^V + ^_£_Ai_^, (4.i7) 

\\2 az A I '^ 2 az A V y ' 

with z being determined by conservation of total mass in the system 

2c + 2z + 4w + 6u = g. (4.18) 

In the case of small grinding, (J3 <C 1), with g and all other parameters being 0(1), we find 

( 2g(3 2 \ 1/3 / g(3 2 \ 1/3 



1 
a 



v 



\3(av + 2 ' \l2(av + 2 ' u lq x 

( H \ 1/3 e (419) 
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In this case most of the mass is in hexamers with a little in tetramers and very little in dimers. 
In the asymptotic limit of a ~ £ 3> 1 and all other parameters 0(1), we find 






1/3 






1/3 



//' 



18£ 



1/3 



Q 

u = -. 

6 



(4.20) 



This differs significantly from the other asymptotic scaling as, not only are c and z both small, 
they are now different orders of magnitude, with c <C z. We next analyse the stability of these 
symmetric states. 



4.2 Stability of symmetric state 



In deriving the above solutions (j4.16p — (I4.17p . we have assumed chiral symmetry, that is, 6 = 
= ip = 0. We now turn to analyse the validity of this assumption. Linearising the system 
of equations fl4.13p - fl4.15p which govern the chiralities, we determine whether the symmetric 
solution is stable from 
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(4.21) 
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For later calculations it is useful to know the determinant of this matrix. Using the steady-state 
solutions (14.161) . the determinant simplifies to 



D 



3c 



{2ac + £z) 2 {a£z 2 -Apfi). 



(4.22) 



For general parameter values, the signs of the real parts of the eigenvalues of the matrix in 
(I4.2ip are not clear. However, using the asymptotic result (I4.19p . for /? < 1, we obtain the 
simpler matrix 



/ 



-P 

'p 2 e{i+av\ 



12 



0V3 (i 



+ au 



1/3 



2/3 



V l2 Q J 
whose characteristic polynomial is 



P 

' P 2 Q(^+av) 
12 

( Pq 2 

\\%{i+av) 



Pi 



\ 



1/3 



1/3 



2 



-W - P 1 '* 



£ + ai> 
2P 2 g 

m+av) 2 ) 



1/3 



2/3 



12 Q 



(4.23) 



= q 3 + fivq 2 + fiv(±p 2 Q(£ + 



1/3 

n//)j q-D, 



(4.24) 



Formally D is the determinant of the matrix in (I4.23p . which is zero, giving a zero eigenvalue, 
which indicates marginal stability. Hence, we return to the more accurate matrix in (I4.2ip . 
which gives D ~ — /3 2 /iz/. The polynomial (I4.24p thus has roots 



Qi 



■fiv, q 2 



'P 2 e {£+au\ 



1/3 



12 



Qs 



' 12/3 4 ' 

A au +0. 



1/3 



(4.25) 
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This means that the symmetric state is always linearly stable for this asymptotic scaling. We 

expect to observe evolution on three distinct timescales, one of 0(1), one of 0((3^ 2 ^ 3 ) and one 

of (9(/3" 4 / 3 ). 

We now consider the other asymptotic limit, namely, a ~ £ ^> 1 and all other parameters 

are 0(1). In this case, taking the leading order terms in each row, the stability matrix in (I4.2ip 

reduces to 

/ «,,,, fi2/3\ 2 / 3 a,.,, fi2/3\ 2 / 3 n \ 



m 

/3 e 2£2\l/3 



<Hf)' 







1/3 



i m) 



(/3e 2 ey /3 



1/3 



■C¥) 

f§ie\ 1/3 



(4.26) 



/ 



which again formally has a zero determinant. The characteristic polynomial is 

= q 3 + q 2 + Qf3fiisq - D, 



(4.27) 



wherein we again take the more accurate determinant obtained from a higher-order expansion 
of (I4.2ip . namely D = /3 2 /iz/. The eigenvalues are then given by 



Qi 



'(3g 2 e\ 1/3 
144 / 



q 2 , 3 ~ ±V/3^ 



'12/3' 



1/3 



(4.28) 



We now observe that there is always one stable and two unstable eigenvalues, so we deduce 
that the system breaks symmetry in the case a ~ £ ^> 1. The first eigenvalue corresponds to a 
faster timescale where t ~ (9(£~ 2 / 3 ) whilst the latter two correspond to the slow timescale where 

t = o(e /3 ). 

4.3 Simulation results 
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Figure 10: Illustration of the evolution of the total concentrations C2,z, w,u for a numerical 
solution of the system truncated at hexamers ( I4.ip -( J4~T1) in the limit a ~ £ ^> 1. Since model 
equations are in nondimensional form, the time units are arbitrary. The parameters are a = 
£ = 30, v = 0.5, (3 = /i = 1, and the initial data is xq(0) = ye(0) = 0.06, £4(0) = 1/4(0) = 0.01, 
£2(0) = 0.051, 1/2(0) = 0.049, c 2 (0) = 0. Note the time axis has a logarithmic scale. 



We briefly review the results of a numerical simulation of f)4.ip -( l4T7|) in the case a ~ £ 3> 1 
to illustrate the symmetry-breaking observed therein. Although the numerical simulation used 
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Figure 11: Graph of the evolution of the chiralities against time on a log- log scale; results of 
numerical simulation of the same hexamer-truncated system, with identical initial data and 
parameters as in Figure [TUJ 

the variables x^ and y^ (k = 2,4,6) and c 2 , we plot the total concentrations z,w,u in Figure 
[TUJ The initial conditions have a slight imbalance in the handedness of small crystals (x 2 ,y 2 ). 
The chiralities of small (x 2 , 2/2, z ), medium (x 4 , 2/4, w), and larger (x 6 , y 6 , u) are plotted in Figure 
[TT1 on a log-log scale. Whilst Figure [JJJ shows the concentrations in the system has equilibrated 
by t = 10, at this stage the chiralities are in a metastable state, that is, a long plateau in the 
chiralities between t = 10 and t = 10 3 where little appears to change. There then follows a 
period of equilibration of chirality on the longer timescale when t ~ 10 4 . We have observed this 
significant delay between the equilibration of concentrations and that of chiralities in a large 
number of simulations. The reason for this difference in timescales is due to the differences in 
the sizes of the eigenvalues in (I4.25P . 

We have also investigated the case (3^1 with all other parameters 0(1) to verify that 
this case does indeed approach the racemic state at large times (that is, 9, <fi, ( — )■ as t — )■ 
00). However, once again the difference in timescales can be observed, with the concentrations 
reaching equilibration on a faster timescale than the chiralities, due to the different magnitudes 
of eigenvalues (14.281) . 



5 New simplifications of the system 



We return to the equations (I2.35l) - fl2.39l) in the case 5 = 0, now writing x 2 = x and y 
obtain 



dc 

dT 
dx 

dt 

f\lj 

— = lie- [ivy - aye + (3(N y - y + y A ) - £y 2 - £yN y , 



dN x 
dt 

dNy 

dt 



— 2/ic + \iv{x + y) — ac(N x + N y ), 

fie — fiux — axe + f3(N x — x + x±) — £x 2 — £xN x , 



fie — \ivx + (3(N X — x) — £,xN x , 
lie - [ivy +P(N y -y)- £yN y , 



y 2 to 

(5.1) 
(5.2) 
(5.3) 
(5.4) 
(5.5) 
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which are not closed, since £4,2/4 appear on the RHS's of (15.21) and (15. 3p . hence we need to find 
formulae to determine x± and 2/4 in terms of x, y, N x , N y . 

One way of achieving this is to expand the system to include other properties of the distri- 
bution of cluster sizes. For example, equations governing the mass of crystals in each chirality 
can be derived as 

-^ = 2fic - 2fiux + 2acN x , -^- = 2\ic - 2fivy + 2acN y . (5.6) 

These introduce no more new new quantities into the macroscopic system of equations, and do 
not rely on knowing 24 or 1/4, (although they do require knowledge of x and y). 

In the remainder of this section we consider various potential formulae for X4, y<± in terms 
of macroscopic quantities so that a macroscopic system can be constructed. We then analyse 
such macroscopic systems in two specific limits to show that predictions relating to symmetry- 
breaking can be made. 

5.1 Reductions 

The equations governing the larger cluster sizes Xk, yk, are 

dx 2k 



dt 



/3(x 2k +2 ~ x 2k ) - (x 2 k ~ x 2k - 2 )(ac + £x); (5.7) 



in general this has solutions of the form x 2k = Y,j Aj{t)Aj~ x : where Aj are parameters (typically 
taking values between unity (corresponding to a steady-state in which mass is being added to 
the distribution) and ac +* x (the equilibrium value); and Aj(t) are time-dependent; for some Aj, 
Aj will be constant. 

We assume that the distribution of each chirality of cluster is given by 

x 2k = xll-y I , I/a* = 3/11 — -^— J , (5.8) 

since solutions of this form may be steady-states of the governing equations (15. 7p . However, in 
our approximations for X4 and 1/4 the parameters X x , X y are permitted to vary with time in some 
way that depends on other quantities in the model equations. The resulting expressions for the 
macroscopic number and mass quantities are 



x / j X 2 k XA X , 
fc=l 


fe=l 


(5.9) 


00 

^2 2kx 2k = 2x\ 2 x , 


00 

Qy = Yl 2k V2k = 2j/Aj. 


(5.10) 



Qx 

fc=l fe=i 

Our aim is to find a simpler expression for the terms X4 and 2/4 which occur in f)5.2l) - fl5.3p . these 
are given by 24 = x(l — 1/X X ) where 

x N* L = _Q^ = JQ£ (5 n) 

x 2N X V2x V ; 

hence 

x 2 2xN x [2x .„ ,„. 

X4 = x — — — , X4 = x , or X4 = x — xJ — . (5-12) 

-''a: Qx V Qx 

There are thus three possible reductions of the equations (j5.ip — (j5.5p . each eliminating one of 
x, N x , g x (and the corresponding y, N y , g y ). We consider each reduction in turn in the following 
subsections. Since some of these reductions involve g x , g y , we also use the evolution equations 
(15. 6p for these quantities. 
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5.2 Reduction 1: to x,y, N x , N y 

Here we assume A^ = N x /x, X y = N y /y, so, in addition to ( 15. ip . (J5.4| - (l5.5p the equations of 
motion are 



jjlc — \ivx + (3N X — — £x — £xN x , (5.13) 



x N x 
dt ^ r "" ' r y N„ s " s "" v "" 



dt ^ n ' x N. 



Hc-/jLvy + /3N y -^-Zy 2 -tyN y ; (5.14) 



we have no need of the densities £> x , g y in this formulation. 

The disadvantage of this reduction is that, due to (15. lip , the total mass is given by 

2N 2 2N 2 
g = 2c+g x + g v = 2c + — - + — y -, (5.15) 

x y 

and there is no guarantee that this will be conserved. 

We once again consider the system in terms of total concentrations and relative chiralities 
by applying the transformation 



x 



§z(l+0), y = \z{l-9), N x = \N{l+<}>), N y = \N{l- 



(5.16) 



to obtain the equations 
dc 



1# —2/j,c + f/,vz — acN, (5-17) 

dz n nnT Pz 2 (l + 8 2 -29(j)) 

- = 2^-^z-acz + PN- N{1 _^ 2) 

-l£z 2 (l + 9 2 )-^zN(l + 94>), (5.18) 

— = 2fic-fiuz + (3N-/3z-^zN{l + 9(f)). (5.19) 

de ( l,*r 2 P Z ld An 

(5.21) 
These equations have the symmetric steady-state given by 9 = = <fi and c, z, N satisfying 

\wz 2(3N(2fi + aN) 



2\x + aN ' (2/3 + £ JV) (2/x + aN) + 2a(ivN ' 



(5.22) 



from (I5.17P and ( I5.19p . Note that the steady state value of N will depend upon the initial 
conditions, it is not determined by f !5.18[) . This is because the steady-state equations obtained 
by setting the time derivatives in (j5.17p ~ (15.19p are not independent. The difference (J5.18p ~ (j5.19p 
is equal to z/N times the sum (J5.17p + (j5.19p . 
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In subsections 15.2.11 and 15.2.21 below, so as to discuss the stability of a solution in the two 
asymptotic regimes /? <C 1 and a ~ £ 3> 1, we augment the steady-state equations (I5.17I) - (I5.19I) 
with the condition g = 2N 2 /z, with g assumed to be 0(1). 

The linear stability of 9 = = </> is given by assuming 9 and <fi are small, yielding the system 




/ (2pc iz_ Bz BN\ (fW Bz _ fiV\ \ 
\ z 2 N z J \ z N 2 ) 



(5.23) 



An instability of the symmetric solution is indicated by the determinant of this matrix being 
negative. Substituting ( 15. 22ft into the determinant, yields 

det = IM^ - ajN 2 ) 

4Bfi + 2aBN + 2^N + 2afiuN + a^N 2 ' v ' ; 



Hence we find that the symmetric (racemic) state is unstable if N > 2w/x/3/a^, that is, large 
aggregation rates (a, £) and slow grinding (B) are preferable for symmetry-breaking. 

We consider two specific asymptotic limits of parameter values so as to derive specific results 
for steady-states and conditions on stability. In both limits, we have that the aggregation 
rates dominate fragmentation (a ~ £ ^> B), so that the system is strongly biased towards the 
formation of crystals and the dimer concentrations are small. In the first case we assume that 
the fragmentation is small and the aggregation rates are of a similar scale to the interconversion 
of dimers (/3^C/i~a~£ = (9(l)); whilst the second has a fragmentation rate of similar size 
to the dimer conversion rates and larger aggregation rates (ct~£;^> ; u~/3 = 0(1)). 

5.2.1 Asymptotic limit 1: f3 <C 1 

In the case of asymptotic limit 1, /3 <C 1, we find the steady-state solution 



m P° 2/3 Bv 

N~a- z ~ , c~- . (5.25) 

From (I5.24p . we find an instability if g > g c := 4/i(^ + ai/)/a^. That is, larger masses (g) favour 
symmetry-breaking, as do larger aggregation rates (a,£). The eigenvalues of (I5.23P in this limit 
are q\ = —\xv - a fast stable mode of the dynamics and 

which indicates a slowly growing instability when g > g c . Hence the balace of achiral to chiral 
morphologies of smaller clusters (u) also influences the propensity for non-racemic solution. 
However, since the dynamics described by this model does not conserve total mass, the results 
from this should be treated with some caution, and we now analyse models which do conserve 
total mass. 

5.2.2 Asymptotic limit 2: a ~ £ > 1 

In this case we find the steady-state solution is given by 

2/3 4/xz/ [J 
z~—, c J — . 5.27 

f a V €e 
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The condition following from H5.24[) then implies that we have an instability if q > q c :— A(a/ql <C 
1. The eigenvalues of the stability matrix are qx = — \\Z(3q£,, which is large and negative, 
indicating attraction to some lower dimensional solution over a relatively fast timescale; the 



eigenvector being (1, 0) T showing that 8 — > 0. The other eigenvalue is q2 = 2(auJ/3/q^ <C 1, and 
corresponds to a slow growth of the chirality of the solution, since it relates to the eigenvector 
(0, 1) T . Assuming the system is initiated near its symmetric solution (# = = 0), this shows that 
the distribution of clusters changes its chirality first, whilst the dimer concentrations remain, at 
least to leading order, racemic. We expect that at a later stage the chirality of the dimers too 
will become nonzero. 



5.3 Reduction 2: to x, y, g x , g y 

Here we eliminate x 4 = x(l — 1/A X .), y 4 = y{\ — 1/A y ) together with N x and N y using 

[XQ~a 



A,. 



Qx » 
2? Xy 



V 



N x 



2 ' y 



VQy 

2 : 



(5.28) 



leaving a system of equations for (c, x, y, g x , g y 
dc 



dt 

dx 
dt 

dg x 
dt 



lau(x + y) - 2/j,c - \f2ac L/xg^ + ^yg~ y 

s- 2 /= / Q% n Q% n 

(ac — (ivx — acx — £x — £xJ — — \~ P\l ~T) P x 




-2(avx + 2 (ac + 2ac 



xg x 



(5.29) 

(5.30) 
(5.31) 



with similar equations for y, g y . Transforming to total concentrations and relative chiralities by 

way of 



x = h(l + 6), y=\z{\ 



g x = 1 1 R(l + 0, g y = lR(l-C), 



we find 



dc 

dt 
dz 

dt 



(AVZ — 2(AC ■ 



acy/zR 

2^ 



^i+e)(i+o + ^(i-6)(i-0 



(5.32) 



(5.33) 



2(xc — (avz — acz — \£ > z 2 (l + 6 2 ) 



/3VzR 
" 2v/2 
Zz^R 1 ' 
4^2 

/3z 3 / 2 



'2R 



^(1+^(1+0 + ^(1-^(1-0 

(l + ^)3/2 (1 + C) l/2 +(1 _^3/2 (1 _ C) lA 
\3/2" 



[l+9f 2 (1 



(l + C )l/2 (1 _ C) l/2 



cLR 
dt" 



2(avz + A(ac + \acy[2zR ^(1 + ^(1+0 + ^(1-^(1-0 



(5.34) 



(5.35) 
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together with the equations fj5.38[) — ( T5.39J) for the relative chiralities 6 and (, which will be 
analysed later. 

Since the equations for dR/ddt and dc/dt are essentially the same, we obtain a third piece 
of information from the requirement that the total mass in the system is unchanged from the 
initial data, hence the new middle equation above. Solving these we find c = \{g — R) and use 
this in place of the equation for c. 

In the symmetric case (8 = ( = 0) we obtain the steady-state conditions 



= 2/j.vz — 4/xc — acv2zR, g = R + 2c, 

= 2uc — uuz — acz — \^z 2 + \(3\/2zR — Bz\ \ — . 

2 z V R 2 V 2 



For small 6, (, the equations for the chiralities can be approximated by 



d0 
dt 



2fjLC 



+ &* + hP\ 



2'- ' 2 \l 2 Z 2 1 V D ' O 




: S 2^2z~R 4V 2 K ' 



dt 



2uvz zR 

R V 2 



/2fiuz 4/ic\ 



(5.36) 



(5.37) 



(5.38) 
(5.39) 



We analyse the stability of the symmetric (racemic) state in the two limits /3 <C 1 and a ~ £ ^> 1 
in the next subsections. 

5.3.1 Asymptotic limit 1: (3 < 1 

In this case, solving the conditions ( 15.36l) -( l5.37l) asymptotically, we find 

2/3 0v 



i + au' 



£ + az/ 



r ~ e - 2c - 



(5.40) 



Substituting these values into the differential equations which determine the stability of the 
racemic state leads to 




/ 



—\xv 



av 



06 \ 



4 V £ + ai/ 
«z//3 3 / 2 



V e(£ + cw) (e + «^) 3/2 v^/ 




(5.41) 



Formally this matrix has eigenvalues of zero and — \xv. Since the zero eigenvalue indicates 
marginal stability of the racemic solution, we need to consider higher-order terms to obtain a 
more definite result. 

Going to higher order, gives the determinant of the resulting matrix as — a^v/{av-\-^) 2 hence 
the eigenvalues are 

9i = -V-v, and 92 = —, —?rz } (5.42) 
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the former indicating a rapid decay of 9 (corresponding to the eigenvector (1, 0) T ), and the latter 
showing a slow divergence from the racemic state in the (^-direction, at leading order, according 
to 

Hence in the case /3 -C 1, we find an instability of the symmetric solution for all other parameter 
values. 

5.3.2 Asymptotic limit 2: a ~ £ > 1 

In this case, solving the conditions ( I5.36p -f l5.37p asymptotically, we find 

z~—, c~ W— , R~g-2c. (5.44) 

4 « V 0? 

Substituting these values into the differential equations f l5.38p -f l5.39p which determine the sta- 
bility of the racemic state leads to 




4/3/iz/ 4/3/xz/ 



\ q£ et I 



(5-45) 



hence the eigenvalues are q\ = — \\/f3g£, and q<i = 4/iz//3/ g£, (in the above o(\/£) means a quantity 
q satisfying q <C \/£ as £ — >■ oo). Whilst the former indicates the existence of a stable manifold 
(with a fast rate of attraction), the latter shows that there is also an unstable manifold. Although 
the timescale associated with this is much slower, it shows that the symmetric (racemic) state 
is unstable. 

5.4 Reduction 3: to N x , N y , g x , g y 

In this case our aim is to retain only information on the number and typical size of crystal 
distribution, so we eliminate the dimer concentrations x, y, using 

A =^- A =-*L x= 2 -^l v= 2 -^l (546) 

X 2iV Xy 2AV g x ' V g y ■ {bAb) 

These transformations reformulate the governing equations fl5.ip - fl5.6p to 

H N N 2 If N 3 

5_!i = 1^-^)+^ -2(^1/ + /?)-£ --A^L, (5.47) 

at g x g x 

dN , N 2 2£W 3 

^ = I /i ( e -J2)+/W y -2fcu/ + ) 8)-*-- L -*-, (5.48) 

Qt gy gy 

^ = (e-R^ + aN,)-^^., (5.49) 

^ = (g-R)^ + aN y )-^l, (5.50) 

where i? := g x + g y . We now transform to total concentrations (N, R) and relative chiralities 
(</> and £) via 

N x =±N(l + <j>), N y = \N{\-4>), q x = \R[1 + Q, g y = \R{l-C), (5.51) 

26 



together with c = \{q — R), to obtain 



dR 



[g-R)(2fj, + aN)- 



4/iz/iV 2 (i + 2 -2<K) 



dN 

— = v( 8 -R)+/3N 
at 

AT 2 r 

2(^+/3)(l + 2 -20C) + ^iV(l + 30 2 -3</)C-0 3 C) 



~dt 



dt 



R(i-Cj 
1 dN 



N dt 

N 



R(l-C 



2(/?+H(20-C-0 2 C) + eiV(30-C + 3 -30 2 O 



a(£-i?)iV0_ ldi? 4/iWV 2 (20-C-0 2 C) 
S fl IP i? 2 (l - C 2 ) 



(5.52) 
(5.53) 

(5.54) 

(5.55) 



We now analyse this system in more detail, since this set of equations conserves mass, and is 
easier to analyse than fl5.33p - fl5.35p due to the absence of square roots. We consider the two 
asymptotic limits (/5<1 and a ~ £ ^> 1) in which, at steady-state, the majority of mass is in 
the form of clusters. 

5.4.1 The symmetric steady-state 

Putting ( = = (f), we find the symmetric steady-state is given by 

4/ii/iV 2 



= 

= 

the former is solved by one of 



[q - R)(2/jl + aN) - 



R 



N 2 £iV 3 
(jl(q -R)+/3N- 2(jjlv + f 3 )^~^Jf 



(5.56) 

(5.57) 



R=h\i± 



WfiuN 2 



\ (2^ + aN)g 2 r 



N = afl(g - R) | ;| 

8yUZ/ 



32/iV 



\ a 2 R(g-R) 



(5.5* 



(5.59) 



More complete asymptotic solutions will be derived in Sections 15.4.31 and 15.4.41 

5.4.2 Stability of the symmetric state 

We now consider the stability of the symmetric steady-state. For small 0, ( we have 



Rd_ 

Ndi 



\c 



^~ R)R 2/m^iv \ 



. 8fwN\ (g-R)(2fjrHxN)R 

a{g-R) — — J Sfiu 



\ 



N 2 



(5.60) 



and this is unstable if the determinant of this matrix is negative. Now we consider the two 
asymptotic limits in more detail. 
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5.4.3 Asymptotic limit 1: (5 <C 1 



When fragmentation is slow, that is, /? -C 1, at steady-state we have iV = 0(\f]3) and i? = 
g — 0{j5). Balancing terms in ( I5.56l) - (l5.57p we find the same leading order equation twice, 
namely 2vN 2 = (3 g(g — R) . Taking the difference of the two yields an independent equation 
from higher order terms, hence we obtain 



N 



(3g 



i + av' 



R ~ Q — 



2vf3 



(5.61) 



Note that this result implies that the dimer concentrations are small, with c ~ z and c ~ 

Substituting these expressions into those for the stability of the symmetric steady-state (j5.60p . 
we find 

This matrix has one stable eigenvalue (corresponding to (1, 0) T and hence the decay of <fi whilst 
( remains invariant), the unstable eigenvector is (1,4) T , hence we find 





exp 






(5.63) 



If we compare the timescale of this solution to that over which the concentrations N, R vary, we 
find that symmetry-breaking occurs on a slower timescale than the evolution of cluster masses 
and numbers. This is illustrated in the numerical simulation of equations (I5.47l) - fl5.50l) shown 
in Figure [12j More specifically, the time-scale increases with the mass in the system, and with 
the ratio of aggregation to fragmentation rates, (ctz/ + £)//?, and is inversely related to the chiral 
switching rate of small clusters (//i / ). 




■ oeoooooooo 
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10° 



Figure 12: Graph of concentrations N x , N y , g x , g y ,c against time on a logarithmic time for the 
asymptotic limit 1, with initial conditions N x = 0.2 = N y , g x = 0.45, g y = 0.44, other parameters 



given by a = 1 = £ = //, (3 = 
the time units are arbitrary. 



0.01 , g = 8. Since model equations are in nondimensional form, 
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5.4.4 Asymptotic limit 2: a ~ £ > 1 

In this case we retain the assumptions that [i, v — 0(1), however, we now impose /3 = 0(1) and 
a ~ £ 3> 1. For a steady-state, we require the scalings N = 0(l/>/D and g — R = 0(l/£ 3 / 2 ). 
Specifically, solving (I5.56p - fl5.57p we find 

AT [PC II ^ V W& (V>RA\ 

N ~ \Hr, R~ Q \~r, (5.64) 

V ? «£ V s 

hence the dimer concentrations c = 5(0 -#) ~ iV 3 = £>(l/£ 3/2 ) and 2 = 2N 2 /g ~ iV 2 = 0(1/0- 
More precisely, c ~ (2/iu/a)Jp/g^ and z ~ 2/3/£, in contrast with the previous asymptotic 
scaling which gave z ~ iV 2 ). 

To determine the timescales for crystal growth and dissolution, we use f)5.64p to define 

N~n(t)jMl, R^g-^^-M (5.65) 

v ag V ? 

and so rewrite the governing equations (I5.52p - fl5.53p as 

d^ = ^( — TE^J' (5 ' 66) 

dt V £ V « V/W 

Here, the former equation for n(t) corresponds to the slower timescale, with a rate /3, the rate 



of equilibration of r(t) being a J /3q/£. 

The stability of the symmetric state is determined by 



Rd_ ( 0(t) \ / -2VM VM\(<f>\ 

Ndt { at) J " I -^Mig ^ [c r [b ^ } 



This matrix has one large negative eigenvalue (~ —2yfj3g£) and one (smaller) positive eigenvalue 
(~ 4//z/); the former corresponds to (1,0) T hence the decay of 0, whilst the latter corresponds 
to the eigenvector (1,2) T . Hence the system (I5.68P has the solution 

( C ) ~ ° ( 2 ) eXP ( 4/iZ/t \/5) ' (5 ' 69) 

The chiralities evolve on two timescales, the faster being 2/3 corresponding to the stable eigen- 



value of (I5.68P and the slower unstable rate being 4/xz/W/3/£^. This timescale is similar to (I5.63p . 
being dependent on mass and the ratio of aggregation to fragmentation, and inversely propor- 
tional to the chiral switching rate of dimers (//f). 

5.5 The asymmetric steady-state 

Since the symmetric state can be unstable, there must be some other large-time asymmetric 

attractor(s) for the system, which we now aim to find. From ( 15.471) and ( 15.491) . at steady-state, 

we have 

4uz/iV 2 N 2 

2c 2 (2/i + aN x ) = -^— £, fic 2 + f3N x = 2(/zz/ + + £N X )-Z. (5.70) 

Qx Qx 
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Figure 13: Graph of the concentrations N x ,N y , g x , g y ,c against time on a logarithmic time for 
the asymptotic limit 2, with initial conditions N x = 0.2 = N y , g x = 0.45, g y = 0.44, other 
parameters given by a = 10 = £, fi = 1 = \x, v — 0.5, g = 2. Since model equations are in 
nondimensional form, the time units are arbitrary. 



Taking the ratio of these we find a single quadratic equation for N x 



= at-Ni 



C--2 



a/3 — apiv — £/i \ N x + /3/j,, 



(5.71) 



with an identical one for N y . Hence there is the possibility of distinct solutions for N x and N y 
if both roots of (15.71 j) are positive; this occurs if 



c 2 < 



/3/iz/ 



a(3 + £fj, + ol[lv + 2i/a/3^/i 
Given N x (N y ), we then have to solve one of ( 15.701) to find g x (g y ), via 

llxvNl 



(5.72) 



Qi 



c 2 (n + aN x ) 



(5.73) 



and then satisfy the consistency condition that g x + g y + 2c 2 = g. After some algebra, this 
condition reduces to 



\a 2 S,c\{(3— a>C2)(g— 2c 2 ) = /3 2 /i 2 z/ 2 — /3/j,vc2[ot/3 + 2a[xv + 2^/j] 

+Hc\[n(ai>+£ t ) 2 + aj3(au—£)]. 



(5.74) 



Being a cubic, it is not straightforward to write down explicit solutions of this equation, hence 
we once again consider the two asymptotic limits (^<1 and a ~ £ ^> 1). 

5.5.1 Asymptotic limit 1: f3 <C 1 

In this case, c 2 = 0(j3) hence we put c 2 = (3C and the consistency condition (15.741) yields 



0((3 3 ) = (3 2 [i/ - (ai/ + £)C] S 



(5.75) 



hence, to leading order, C = v j[av + £) . Unfortunately, the resulting value for c 2 leads to all 
the leading order terms in the linear equation f !5.71|) for N x to cancel. We thus have to find 
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higher order terms in the expansion for c 2 ; due to the form of (j5.75p . the next correction term 
is C(/3 3 / 2 ). Putting c 2 = f3C(l + Cyflf), we find 

~ 2 _ at, [a£Q + ±ii{av + £)] . . 

In order to satisfy the inequality (15.72]) . we require the negative root, that is, C < 0. 
Although the formulae for N x , N y are lengthy, their sum and products simplify to 



a? y at; 

The chirality <fi can be simplified using <ft 2 = 1 — 4I1/S 2 which implies 

2 _ agj - 4/j(q^ + £) 
a Qi + 4/i(az/ + £) 



N x + N y = ^^^f^l, U = N x N y = ^. (5.77) 



(5.78) 



Hence we require g > g c := 4/z(ai/ + £)/«£ in order for the system to have nonsymmetric steady- 
states, that is, the system undergoes a symmetry-breaking bifurcation as g increases through 
g = g c . As the mass in the system increases further, the chirality approaches (±) unity, 
indicating a state in which one handedness of crystal completely dominates the other. 

5.5.2 Asymptotic limit 2: a ~ £ 3> 1 

In this case, the left-hand side of the consistency condition ( 15.74]) is 0{a 2 ^c^) whilst the right- 
hand side is 0(1) + 0(ac 2 l ), which implies the balance c 2 = 0(£~ 3 / 2 ). Solving for c 2 leads 
to 

c 2 ~ aa-M. (5.79) 

a v gt, 
The leading order equation for N x , N y is then 



= a£N 2 - aNJ\(5gi + /3/i, (5.80) 



hence we find the roots 



Wg 2ii [f3 2// 

i\ x , i\ v ~ t — , — w— — , g x , gv ~ g, — • (5.81) 

Since we have either g x ^> N x ^> g y ^> N y or g y ^ N y ^ g x ^ N x , in this asymptotic limit, 
the system is completely dominated by one species or the other. Putting £ = N x + N y and 
II = N x N y we have 2 = 1 - 4II/E 2 ~ 1 - S/i/ag. 

6 Discussion 

We now try to use the above theory and experimental results of Viedma [29] to estimate the 
relevant timescales for symmetry-breaking in a prebiotic world. Extrapolating the data of time 
against grinding rate in rpm from Figure 2 of Viedma [29] suggests times of 2 x 10 5 hours using 
a straight line fit to log(time) against log(rpm) or 1000-3000 hours if log(time) against rpm or 
time against log(rpm) is fitted. A reduction in the speed of grinding in prebiotic circumstances 
is expected since natural processes such as water waves are much more likely to operate at the 
order of a few seconds" 1 or minutes -1 rather than 600 rpm. 
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Similar extrapolations on the number and mass of balls used to much lower amounts gives a 
further reduction of about 3, using a linear fit to log(time) against mass of balls from Figure 1 
of Viedma [25]. There is an equally good straight line fit to time against log (ball- mass) but it is 
then difficult to know how small a mass of balls would be appropriate in the prebiotic scenario. 
There is an additional factor due to the experiments of Viedma being on a small volume of 10 
ml, whereas a sensible volume for prebiotic chemistry is 1000 1, giving an additional factor of 
10 5 . Combining these three factors (10 3 , 3, and 10 5 ) with the 10 days of the original experiment, 
we estimate that the timescale for prebiotic symmetry breaking is (9(3 x 10 9 ) days, which is 
equivalent to the order of about ten million years. 

This extrapolation ignores the time required to arrive at the initial enantiomeric excesses of 
5% used by Viedma [29J from a small asymmetry caused by either a random fluctuation or by the 
parity- violation. Although the observed chiral structures are the minimum energy configurations 
as predicted by parity violation, there is an evens probability that the observed handedness could 
simply be the result of a random fluctuation which was amplified by the same mechanisms. In 
order to perform an example calculation, we take a random fluctuation of the size predicted by 
parity violation, which is of the order of 10 -17 , as suggested by Kondepudi & Nelson [16] . Our 
goal is now to find the time taken to amplify this to an 0(1) (5%) enantiomeric excess. 

The models derived in this paper, for example in Section 15.4.41 predict that the chiral excess 
grows exponentially in time. Assuming, from (I5.69p . that 0(to) = 10 -17 and 4>(ti) = 0.1, then 
the timescale for the growth of this small perturbation is 

ti — to = 

Since the growth of enantiomeric excess is exponential, it only takes 16 times as long for the 
perturbation to grow from 10~ 17 to 10 _1 as from 10 _1 to 1. Hence we only need to increase our 
estimate of the timescale by one power of ten, to 100 million years. 

This estimate should be taken as a very rough estimate, since it relies on extrapolating results 
by many orders of magnitude. Also, given the vast differences in temperature from the putative 
subzero prebiotic world to a tentative hot hydrothermal vent, there could easily be changes in 
timescale by a factor of several orders of magnitude. 

7 Conclusions 

After summarising the existing models of chiral symmetry-breaking processes we have systemat- 
ically derived a model in which through aggregation and fragmentation chiral clusters compete 
for achiral material. The model is closed, in that there is no input of mass into the system, 
although the form of the aggregation and fragmentation rate coefficients mean that there is an 
input of energy, keeping the system away from equilibrium. Furthermore, there is no direct in- 
teraction of clusters of opposite handedness; rather just through a simple competition for achiral 
substrate, the system can spontaneously undergo chiral symmetry-breaking. This model helps 
explain the experimental results of Viedma [21] and Noorduin et al. |21j . 

The microscopic model originally derived has been simplified successively to a minimalistic 
model, which, numerical results show, exhibits symmetry-breaking. Even after this reduction, 
the model is extremely complex to analyse due to the large number of cluster sizes retained 
in the model. Hence we construct two truncated models, one truncated at tetramers, which 
shows no symmetry-breaking and one at hexamers which shows symmetry-breaking under certain 
conditions on the parameter values. Alternative reductions are proposed: instead of retaining 
the concentrations of just a few cluster sizes, we retain information about the shape of the 
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distribution, such as the number of clusters and the total mass of material in clusters of each 
handedness. These reduced models are as simple to analyse as truncated models yet, since 
they more accurately account for the shape of the size-distribution than a truncated model, are 
expected to give models which more easily fit to experimental data. Of course, other ansatzes 
for the shape of the size distributions could be made, and will lead to modified conditions for 
symmetry-breaking; however, we believe that the qualitative results outlined here will not be 
contradicted by analyses of other macroscopic reductions. 

One noteworthy feature of the results shown herein is that the symmetry-breaking is inher- 
ently a product of the two handednesses competing for achiral material. The symmetry-breaking 
does not rely on critical cluster sizes, which are a common feature of theories of crystallisation, 
or on complicated arguments about surface area to volume ratios to make the symmetric state 
unstable. We do not deny that these aspects of crystallisation are genuine, these features are 
present in the phenomena of crystal growth, but they are not the fundamental cause of chiral 
symmetry-breaking. 

More accurate fitting of the models to experimental data could be acheived if one were to 
fit the generalised Becker-Doring model (I2.lip -f l2.12p with realistic rate coefficients. Questions 
to address include elucidating how the number and size distribution at the start of the grinding 
influences the end state. For example, if one were to start with a few large right-handed crystals 
and many small left-handed crystals, would the system convert to entirely left- or entirely right- 
handed crystals ? Answers to these more complex questions may rely on higher moments of the 
size distributions, surface area to volume ratios and critical cluster nuclei sizes. 
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A General theory for crystallisation and grinding with 
competition between polymorphs 

This model can be generalised so as to be applicable to the case of grinding a system undergoing 
crystallisation in which several polymorphs of crystal nucleate simultaneously. It may then be 
possible to use grinding to suppress the growth of one polymorph and allow a less stable form 
to be expressed. In this case, the growth and fragmentation rates of the two polymorphs will 
differ, we denote the two polymorphs by x and y. In place of a, b, a, £, (3 we have a XjT ., a y>r , b x>r , 
a x , r , etc. Hence in place of fl2.20p -f l2.27p we have 



dx r 

"dT 



&x,r—lCl3Ji — 1 Ox,r%r ^x,rC\X r -rOa^r+l-^r+l Px,r%r ~~i~ Px,r+2%r+2 

+ (a Xj r-2C2 + £,x,r-2X2)Xr-2-(a x ,rC2 + £,x,rX2)Xr, (r>4), (Al) 

a y,r-lC\y r -l ~ Oy,rVr ~ Q>y,r c lVr + 0y,r+l3/r+l ~~ Py,rVr + Py,r+2Vr+2 

+ {ay,r-2C2 + €y, r -2y2)y r -2-{0ly,rC2 + £y,ry2)yn ( r > 4 )? ( A2 ) 
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dx 



2 



dt 



dy 2 
dt 



HxC2 ~ [ l x v x^2 — a x.2,C\X 2 + ^,3^3 — («:r,rC 2 + £x,rX2)X r 

oo oo 

+/3 Xj4 X 4 + J2 Px,rX r ~ Yl ^x,kX 2 X k , (A3) 

fc=4 fc=2 

/i y c 2 - H y v y y 2 - a y ^c 1 y 2 + b y , 3 y 3 - (a Vyr c 2 + iy^Vr 

oo oo 

+I3 y ,m + Yl PvrVr - J2 iyMViy^ (A4) 

fc=4 fc=2 



dx 



:-! 



dt 



ax,2^2Cl-fe x ,3^3-a a: ,3ClX3 + 6 a;i 4a;4-(ax,3C2 + ^,3^2)^3 + /3 ;C ,5^5, 



(A5) 



-jT = «y,2l/2Cl-6j / ,3Z/3-S,3ClZ/3 + &|/,4Z/4-(a^3C2+^,3l/2)l/3 + /3y,5?/5, 



dt 

dc 2 

"dT 

dci 



(A6) 

oo 
VxV x X 2 +IJ, y V y y2-(lJ lx +IJ / y)c2^cl^C2-Y C 2( a x,r X r J ray,ryr), (A7) 



fc=2 



2ec 2 - 2<fci - ^ (ctx.fcCiXfc - fc^+i^fc+i + a y ^c x y k - by,k+iyk+i ] 

fc=2 



(A8) 



For simplicity let us consider an example in which all the growth and fragmentation rate 
parameters are independent of cluster size, (a x>r = a x , £ yjt . = £ y , etc. for all r). The thermody- 
namic stability of the two types of crystal depends on their relative interactions with monomers 
from solution, that is, if a x /b x > a y /b y then X is the more stable form. This is because, in the 
absence of c 2 , we can define free energy functions 



a, ^ r ~ 1 




r-l 



Qr = [f) , <#=(?) , (A9) 



which generate the equilibrium distributions 



„eqx _ n xr _ _^£ / Q x c l V eqy _ n y r _ \_ I a y C l 



r 



cr= Q* c r "JL rp > c eiy = Q yr ^ rp (AlQ) 

' a x V b x J a y \ b y J 

If a x /b x < dyjby then the latter (Y) will be the dominant crystal type at equilibrium, whilst X is 
the less stable morphology at equilibrium. These last two words are vital, since, at early times, 
the growth rates depend on the relative sizes of the growth rates a x and a y . It is possible for the 
less stable form to grow first and more quickly from solution, and be observed for a significant 
period of time, since the rate of convergence to equilibrium also depends on the fragmentation 
rates and so can be extremely slow (see Wattis [30] for details). 

In the presence of grinding, the crystal size distributions also depend upon the strength 
of dimer interactions, that is, the growth rates ot x c 2 + £, x x 2l a y c 2 + ^ y y2 and the grinding rates 
/3 X , fly. The steady-state size distributions will depend on the relative growth ratios due to 
grinding (a x c 2 + ^, x x 2 )/(3 x and (a y c 2 + Cyl/2)/Py as well as the more traditional terms due to 
growth from solution, namely a x c\/b x and a y ci/b y . Such systems with dimer interactions have 
been analysed previously by Bolton & Wattis [3]. The presence of dimer interactions can alter the 
size distribution, and in non-symmetric systems such as those analysed here, dimer interactions 
can alter the two distributions differently. Two points are worth noting here: 
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(i) for certain parameter values, the less stable stable form (Y, say, with a y /b y < a x /b x ) may 
be promoted to the more stable morphology by grinding (if {ot y C2+£ y y2) I 'P y is sufficiently greater 
than (a x c 2 +£ x x 2 )/[3 x ); 

(ii) grinding may make a less rapidly nucleating and growing form (Y, say, with a y < a x ) into 
a more rapidly growing form if a y C2+£, y y2 is sufficiently greater than a x C2+£,2X2- 

In systems which can crystallise into three or more forms, we may have the case where x 
is more stable than y and y is more stable than z; thus, at equilibrium x will be observed. 
Furthermore, if a x < a y > a z we may observe type y at early times due to it having faster 
nucleation and growth rates than x and z. However, it is possible that the presence of grinding 
could suppress both x and y and allow z to be expressed, if some combination of the inequalities 

a z C2 + i z Z2 a y c 2 + i y y2 a x c 2 + £ x x 2 , . , 

~T Z > ~T y ' K ' (A11) 

a z > a y , a x , £ 2 > £ x , £ y hold. 
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